Ultracold atoms confined in an optical lattice plus parabolic potential: a closed-form approach 
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We discuss interacting and non-interacting one dimensional atomic systems trapped in an optical lattice plus 
a parabolic potential. We show that, in the tight-binding approximation, the non-interacting problem is exactly 
solvable in terms of Mathieu functions. We use the analytic solutions to study the collective oscillations of 
ideal bosonic and fermionic ensembles induced by small displacements of the parabolic potential. We treat the 
interacting boson problem by numerical diagonalization of the Bose-Hubbard Hamiltonian. From analysis of 
the dependence upon lattice depth of the low-energy excitation spectrum of the interacting system, we consider 
the problems of "fermionization" of a Bose gas, and the superfiuid-Mott insulator transition. The spectrum of 
the noninteracting system turns out to provide a useful guide to understanding the collective oscillations of the 
interacting system, throughout a large and experimentally relevant parameter regime. 
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I. INTRODUCTION 

In recent experiments , quasi-one dimensional 

systems have been realized by tight confinement of gases in 
two dimensions. Due to the enhanced importance of quan- 
tum correlations as dimensionality is reduced, such systems 
exhibit physical phenomena not present in higher dimensions, 
such as "fermionization" of a Bose-Einstein gas. 

The degree of correlation is measured by the ratio of the 
interaction energy to the kinetic energy, 7. For 7 « 1 the 
system is weakly interacting, and at zero temperature most 
atoms are Bose-condensed. In this limit quantum correlations 
are negligible and the dynamics is governed by the mean-field 
Gross-Pitaevski equation. For 7 » 1, the system is highly 
correlated and becomes fermionized, in that repulsive inter- 
actions mimic the effects of the Pauli exclusion principle. In 
this "Tonks-Girardeau" regime, 0, 0], the low energy exci- 
tation spectrum of a bosonic gas resembles that of a gas of 
non-interacting fermions. 

To date, one dimensional systems have been obtained by 
loading a Bose-Einstein condensate into a two-dimensional 
optical lattice which is deep enough to restrict the dynamics 
to one dimension. This procedure creates an array of indepen- 
dent ID tubes. In most experiments, a weak quadratic poten- 
tial is superimposed upon the lattice, in order to confine the 
atoms during the loading proccess. The combined presence 
of the periodic and quadratic potentials substantially modi- 
fies the dynamics of the trapped atoms compared to the cases 
when only one of the two potentials is presen t as shown both 
experimentally and theoretically fl8ll9t ll0l.ll ill 1211 1311 . 

In this paper we study both ideal and interacting bosonic 
systems in such potentials. In contrast to previous stud- 
ies of ideal systems, which used numerical 1 10] or approxi- 
mate solutions |H[nJ[2>Q> nere we show that the single- 
particle problem is exactly solvable in terms of Mathieu 
functions 1 14 115111611 . We use analytic solutions to fully char- 
acterize the energy spectrum and eigenfunctions in the vari- 
ous regimes of the trapping potentials and to provide analytic 
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expressions for the oscillations of the center of mass of both 
ideal bosons and fermions that are induced by small displace- 
ments of the parabolic trap. These expressions for the dipolar 
motion may be tested in experiments. 

We further analyze the low-energy spectrum of interact- 
ing bosons by means of exact diagonalizations of the Bose- 
Hubbard Hamiltonian, and identify the conditions required 
for fermionization to occur. Moreover, we show that spe- 
cific changes in the spectrum of the fermionized system can 
be used to describe the characteristics of a Mott insulator state 
with unit filling at the trap center. 

Center-of-mass oscillations are also studied in the weakly 
interacting and fermionized regimes by comparing exact 
numerical solutions for the interacting system to solutions for 
ideal bosons and fermions, respectively. Because fermioniza- 
tion occurs over a large range of trap parameters, knowledge 
of the properties of the single-particle solutions turns out to 
provide useful insights in the understanding of the complex 
many-body dynamics. In addition, we numerically analyze 
the distribution of frequencies pertaining the modes excited 
during the collective dynamics, for different values of 7. 
This helps us gain qualitative insight in the dynamics even in 
the intermediate regime where interaction and kinetic ener- 
gies are comparable and no mapping to ideal gases is possible. 

The presentation of the results is organized as follows: 
In Sec. II we discuss the analytic solutions that describe the 
single-particle physics. These show the existence of two types 
of behavior, according to whether the energy is dominated 
by site hopping or parabolic contributions. Different asymp- 
totic expansions of the eigenfunctions and eigenvalues apply 
to these two regimes, and can be effectively combined to de- 
scribe the full spectrum. 

In Sec. Ill, we then apply the analytic solutions and asymp- 
totic expansions to the description of the collective dynamics 
of non-interacting ensembles of bosons and fermions subject 
to a sudden displacement of the parabolic trap. In contrast 
to the well-known case of the displaced harmonic oscillator, 
a non-trivial modulation of the center of mass motion occurs 
due to the presence of the lattice potential. We derive explicit 
expressions for the initial decay of the amplitude of the oscil- 
lation, to which we refer as effective damping. 

In Sec. IV A, the low-energy spectrum of interacting bosons 
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is studied as a function of the lattice depth. In particular, we 
follow the evolution of the spectrum from ideal bosonic to 
ideal fermionic as the lattice deepens, by means of numeri- 
cal diagonalization of the Hamiltonian for a moderate number 
of atoms and wells. We specify the necessary conditions for 
the formation of a Mott state at the center of the trap, and use 
the Fermi-Bose mapping to link its appearance and the reduc- 
tion of number fluctuations to the population of high-energy 
localized states at the Fermi level. 

Section IV B is dedicated to the study of the collective 
dipole dynamics of an interacting bosonic gas by numerical 
calculations of the exact quantal dynamics. The dynamics is 
studied for two different scenarios that are experimentally re- 
alizable. The first corresponds to experiments in which the 
trapping potentials are kept fixed and the interatomic scatter- 
ing length is varied (e.g. by use of a Feshbach resonance), 
Sec. IV Bl. Here, parameters have been chosen such that no 
Mott insulator at the center of the trap is formed in the large 
7 limit. The second scenario corresponds to experiments in 
which the lattice depth is increased while the frequency of the 
parabolic potential is kept fixed, Sec. IV B2. In this case a unit 
filled Mott insulator is formed at the trap center and the inhi- 
bition of the transport properties of the system is observed as 
the lattice deepens. 



binding approximation yields the following equations of mo- 
tion for the amplitudes {zj(t)}: 
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where J is the tunneling matrix element between nearest 
neighboring lattice sites. In Eq.(|3} the overall energy shift 
e given by 



dxwQ(x)H wo(x)dx, 



(7) 



has been set to zero. 



II. SINGLE PARTICLE PROBLEM 

A. Tight binding solution 

The dynamics of a single atom in a one dimensional op- 
tical lattice plus a parabolic potential is described by the 
Schrodinger equation 
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where <£>(a;) is the atomic wave function, lot is the trapping 
frequency of the external quadratic potential, V a is the optical 
lattice depth which is determined by the intensity of the laser 
beams, a = A/2 is the lattice spacing, A is the wavelength of 
the lasers and m is the atomic mass. 

If the atom is loaded into the lowest vibrational state of 
each lattice well and the dynamics induced by external per- 
turbations does not generate interband transitions, the wave 
function $(a;) can be expanded in terms of first-band Wannier 
functions only iflTi fTsIl 



(t)w (x- ja), 



(2) 



where wo(x — ja) is the first-band Wannier function centered 
at lattice site j, t is time, and {zj (t)} are complex amplitudes. 
For a deep enough lattice, tunneling to second nearest neigh- 
bors can be ignored. This approximation known as the tight 



B. Stationary solutions 



The stationary solutions of Eq.(|3j are of the form (t) = 

fj > e~ lE " t ^ h , with E„ and /j n ' the n th eigenenergy and 
eigenstate, respectively. Substitution into Eq.(|3} yields 



E n f. 



(n) _ 
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Equation is formally equivalent to the recursion relation 
satisfied by the Fourier coefficients of the periodic Mathieu 
functions with period tt. Therefore the eigenvalue problem 
can be exactly solved by identifying the amplitudes fj' and 
eigenenergies E n with the Fourier coefficients and character- 
istic values of such functions , respectively 0. In terms of 
Mathieu parameters the symmetric and antisymmetric solu- 
tions are 



Jn=2r) 
J 3 

(n=2r+l) 



2tt 



ce2r (x, —q) cos(2jx)dx, (9) 



se 2r (x, —q) sm(2jx)dx, (10) 



E n =2r — 'J a 2r (?) S n= 2 r +1 = "J&2r (?) j (1 1) 



with r — 0,1,2... and cei r {x,q) and se2 r (x,q) the even 
and odd period ir solutions of the Mathieu equation with 
parameter q = j[ and characteristic parameter tt2 r (?) and 

b 2r (q) respectively: d2 ^i r + (a 2r {q) - 2q cos(2a;))ce2 r = 

and ^j^- + (b 2r (q) - 2qcos(2x))se 2r = 0. 
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Solutions of the single particle problem are entirely 
determined by the parameter 9 = 4r' wn i cn i s proportional 
to the ratio of the nearest-neighbor hopping energy J to the 
energy cost f2 for moving a particle from the central site to 
its nearest-neighbor. The eigenvalues and eigenstates of the 
system are in general complicated functions of q. However, 
asymptotic expansions exist in literature which can help 
unveil the underlying physics for different values of q. Most 
of the asymptotic expansions have been available for almost 
50 years due to the work of Meixner and Schafke Hal . In 
the remainder of this section we introduce such asymptotic 
expansions and use them to describe the physics of the system 
in the tunneling dominated regime where q > 1, which we 
call high q regime, and the q < 1, or low q regime. 



1. High q regime (4 J > fi) 



are localized on the sides of the potential. The existence of 
localized and extended states has been tested exp erimentally 
H H and studied theoretically O [H d [Ij. Here we 
show how the asymptotic expansions of the Mathieu solutions 
can be used to characterize them quantitatively. 

Low-energy modes (n <C y/q) in the high q regime 
In the LE regime the average hopping energy J is larger 
than fl. In this regime the eigenmodes have been shown to 
be approximately harmonic oscillator eigenstates 
1 1 311 . Asymptotic expansions valid to describe LE eigenmodes 
in the high q limit have been studied in detail in Ref. (k|, 
where the eigenstates are described in terms of generalized 
Hermite polynomials. Here we simply outline the basic results 
and refer the interested reader to 111 611 and references therein 
for further details. 

The asymptotic expansions for the eigenmodes can be writ- 
ten as 



Most experiments have been developed in the parameter 
regime where q ^> 1. For example for 87 Rb atoms trapped 
in a lattice with A = 810 nm, a value of q > 10 is obtained for 
a lattice depth of 2 Er if cot < 2tt x 538 Hz and for a lattice 
of 50 Efl if lot < 2ir x 7.5 Hz. Here Er is the photon re- 
coil energy Er = h 2 /(2m\ 2 ), corresponding to a frequency 
Er/H — 3.47 kHz. Throughout this paper, in our examples 
we use atoms with the mass of 87 Rb. 

In the high q limit the periodic plus harmonic potential 
possesses two different classes of eigenstates depending 
on their energy. In fact, as we will show, eigenmodes 
can be classified as low or high-energy depending on the 
quantum number n being smaller or larger than 2 1 1 \Jq/2\ | , 
respectively, where ||x|| denotes the closest integer to x. 
Physically, this classification depends on which one of the 
two energy scales in the system, the tunneling or the trapping 
energy, is dominant. To the energy classification corresponds 
a classification based on localization of the modes in the 
potentials. In particular, the low-energy (LE) modes are 
extended around the trap center and high-energy (HE) modes 
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and A n a normalization constant. Notice 
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that the coefficients h, and h, are related to the Hermite 

( r )^2fc 



(2fc+l)!(i--fc)! 



polynomial H n (x) by the relations H2 r (x) = J2k=o x 
and H 2r +i{x) = J2k=o^ 
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The eigenenergies are approximately given by 



J 



e: 



4 



-r\-2q + 4^(r 
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2' 8 2 7 ^/q 



oil 
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If one neglects corrections of order 1/ \fq and higher in Eqs. 
dl 2t and (1131 and keeps only the first two terms in Eq. dl4t . 
the expressions for the eigenmodes and eigenenergies reduce 
to 

En - -n q /2 + n^(n + l-) 



f 



(n) 



V2 



■. exp 



I! 

2 



(15) 



(16) 



V 2 n n\ yqir 2 

The above expressions correspond to the eigenvalues and 



eigenenergies (shifted by —Ctq/2) of a harmonic oscillator 
with an effective trapping frequency u* and an effective mass 
m*. The effective frequency and mass are given by 



hw* 



rn 



fl^/q = hwT 

h 2 

2Ja 2 ' 



(17) 



(18) 



The harmonic oscillator character of the lowest energy modes 
in the combined lattice plus harmonic potential is consistent 
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with the fact that near the bottom of the Bloch band the disper- 
sion relation has the usual free particle form with m replaced 
by m* . It is important to emphasize that the expressions for 
the effective mass and frequency Eqs. d!7i and Jl 8I > are only 
valid in the tight-binding approximation. 

Higher order terms introduce corrections to these harmonic 
oscillator expressions, that become more and more important 
as the quantum number n increases. These corrections come 
from the discrete character of the tight-binding equation. They 
can be calculated by replacing /j™" 1 with f^ n \£,) and taking 
the continuous limit of the hopping term proportional to J in 
Eq.(|8j. This procedure yields: 

(-\^L I + ... + 1? 2 ) /■(») - e f {n \m 

\ 2de i2a ho *de + 2 l; ) s n/ ( y) 

where £ = jfi ee £ , a ho = = (q/i)^ is 
I 



with A n normalization constants. We note that the asymp- 
totic expansions Eas. (l20> and (12 1 i for the HE-modes in the 
high q regime are identical to the asymptotic expansions for 
the modes in the low q regime reported below. While the lat- 
ter are well known in literature, to our knowledge they have 
never been used to describe the high q regime. We actually 
found that as long as n ^fq these asymptotic expansions re- 
produce reasonably well the exact results, even for very large 
q. An example of this is given in Figs. ^and|2] which are 
discussed at the end of this section. 

The high energy eigenstates are almost two-fold degenerate 
with energy spacing mostly determined by Q. In Ref. lll2ll 
the authors show that the localization of these modes can be 
understood by means of a simple semiclassical analysis. By 
utilizing a WKB approximation, the localization of the modes 
can be linked to the appearance of new turning points in 
addition to the classical harmonic oscillator ones for energies 
greater than 2 J. While classical harmonic oscillator turning 
points are reached at zero quasimomentum, the new turning 
points appear when the quasimomentum reaches the end of 
the Brillouin zone and can therefore be associated with Bragg 
scattering induced by the lattice. 

Intermediate states (n ~ ^fq) in the high q regime 



a characteristic length of the system in lattice units, which 
can be understood as an effective harmonic oscillator length, 
and E n = {E n + 2J)/{hw*). To zeroth order in l/y/q, 
iflho —* 00 X the differential equation dl9l reduces to the 
harmonic oscillator Schrodinger equation. Higher order 
corrections given in Eas.( ll2> . d 1 3I > and dl4> . can be calculated 
by treating the higher order derivatives in Eq.(ll9> as a 
perturbation. 

High-energy modes ( n 3> y/q) in the high q regime 

High energy modes are close to position eigenstates since 
for these states the kinetic energy required for an atom to hop 
from one site to the next becomes insufficient to overcome 
the potential energy cost By using asymp- 

totic expansions of the characteristic Mathieu values and func- 
tions, we obtain the following expressions for the spectrum 
and eigenmodes 



(20) 



(21) 

I 

In order to reproduce accurately the energy spectrum in the 
intermediate regime one may expect that many terms in the 
asymptotic expansions have to be kept. To estimate the en- 
ergy range where the spectrum changes character from low 
to high, we solve for the smallest quantum number n c whose 
energy calculated by using the low energy asymptotic expan- 
sion is higher than the one evaluated by using the high energy 
expansion. That is 

E l ™_i > E^ h , (22) 

where n c is required to be even. 
Solution of Eq. d22l gives 

n c w2||Vg72||, (23) 

where ||a;|| denotes the closest integer to x. 

By comparison with the numerically obtained eigenvalues, 
we actually find that using Eq. MAI for n < n c — 1 and 
Eq. (I20> for n > n c — 1 is enough to reproduce the entire 
spectrum quite accurately. Moreover, since at n c the energy is 
approximately given by E nc w 2 J, our analytical findings for 
the transition between harmonic oscillator-like and localized 
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eigenstates are in agreement with the approximate solutions 
found in 1 1 it fl2fl by using a WKB analysis. 
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FIG. 1: Upper panel (a): Spectrum of a particle in combined 
quadratic and periodic potentials as a function of the quantum num- 
ber n. The quadratic trap frequency and the depth of the periodic 
potential are lot — 2tt x 60 Hz and V = 7AEr, respectively. 
Points are numerically obtained values, while crosses are asymptotic 
expansions of the Mathieu characteristic parameters. The arrow in- 
dicates the critical value n c w 2||y q/2\\. Lower panel (b): Energy 
difference DE n between the numerically obtained eigenvalues and 
the asymptotic expansions. 



FIG. 2: Eigenmodes of a particle in a combined quadratic and peri- 
odic potentials as a function of the lattice site j. The quadratic trap 
frequency and the depth of the periodic potential are lot = 2tt x 60 
Hz and V = 7AEr, respectively. Triangles and dots are numeri- 
cally obtained values ("ex" in the legend), while boxes and crosses 
are asymptotic expansions of the Fourier coefficients of the periodic 
Mathieu functions ("as" in the legend). In the upper panel, triangles 
and boxes refer to the n — 2 mode, while dots and crosses refer to 
n = 5. In the lower panel, triangles and boxes refer to the n = 34 
mode, while dots and crosses refer to n = 42. 



In Figs. H an d|2l we compare the above asymptotic approx- 
imations to exact numerical results for the eigenenergies and 
eigenfunctions. The parameters used for the plots are those 
of a system with a lattice depth of IAEr and quadratic trap 
frequency lot — 2tt x 60Hz. These values correspond to 
J = 0M57E R , n = 0.0009^ and q = 157. Figure [T] 
upper panel, shows the lowest 35 eigenenergies as a function 
of the quantum number n. The value n c , which is equal to 
16 in this case, is indicated by an arrow. The crosses rep- 
resent the asymptotic solutions, Eq. 1141 and i20\ . and the 
dots the numerically obtained eigenvalues. On the scale of 
the graph there is essentially no appreciable difference be- 
tween the two solutions for the entire spectrum. The differ- 
ence between the the numerically obtained energies and the 
asymptotic expansions is plotted in the lower panel of Fig. [2 
In the upper panel of Fig|2 asymptotic expressions for the 
LE eigenvectors n = 2(boxes) and n = 5(crosses) are com- 
pared to the numerically obtained eigenmodes (triangles and 
dots respectively). The modes clearly exhibit an harmonic 
oscillator character, and the agreement between the asymp- 
totic and numerical solutions is very good. In the lower panel, 
the n = 34 and n — 42 eigenstates belonging to the region 



n > n c are depicted. These states are localized far from the 
trap center. While the overall shape of the modes is well repro- 
duced by the asymptotic solutions, for the chosen values of n 
small differences between the asymptotic (boxes and crosses 
respectively) and numerical solutions (triangles and dots re- 
spectively) can be observed. As expected, the convergence 
of the asymptotic expansion to the exact solution is better for 
n = 42 than for n — 34, as the former has a larger value of n 
than the latter. 



2. Low q regime (4 J < Q) 

This parameter regime is relevant for deep lattices. When 
4 J < fl the kinetic energy required for an atom to hop from 
one site to the next one is insufficient to overcome the trapping 
energy even at the trap center and all the modes are localized. 
This is consistent with the previous analysis, because when 
4 J < O, n c is less than one. The asymptotic expressions that 
describe this regime are 11411 
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FIG. 3: Upper panel (a): Spectrum of a particle in combined 
quadratic and periodic potentials as a function of the quantum num- 
ber n. The quadratic trap frequency and the depth of the periodic po- 
tential are ujt — 2n x 60 Hz and V = 50.0Er, respectively. Points 
are numerically obtained values, while crosses are asymptotic ex- 
pansions of the Mathieu characteristic parameters in the low q limit. 
Lower panel (b): Energy difference DE„ between the numerically 
obtained eigenvalues and the asymptotic expansions. 
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with A„ normalization constants. In Fig. [3]the above expan- 
sions for the energies are compared with the numerically cal- 
culated spectrum. Here the lattice is 50-Er deep and the exter- 
nal trap frequency is lot = 2tt x QOHz. These lead to values 
of J, il and q given by J = 2.9 x Kr 5 .E fl ,, D, = 0.0009£: fl 
and q = 0.13. The asymptotic and numerical solutions per- 
fectly agree on the scale of the graph. The energy difference 
between the numerically obtained eigenvalues and the asymp- 
totic expansions is plotted in Fig. [3] lower panel. 



III. CENTER OF MASS EVOLUTION OF A DISPLACED 
SYSTEM 
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In recent experiments, the transport properties of one di- 
mensional Bose-Einstein condensates loaded in an optical lat- 
tice have been studied after a sudden displacement of the 
quadratic trap |19[. A strong dissipative dynamics was ob- 
served even for very small displacements and shallow depths 
of the optical lattice. This should be contrasted with previ- 
ous experiments performed with weakly interacting 3D gases 
where very small damping of the center of mass motion was 
observed for small trap displacements 12(1 Elll . Recent the- 
oretical studies have demonstrated that the strongly damped 
oscillations observed in one dimensional systems reflect the 
importance of q uantum fluctuations as the dimensionality is 
reduced E3 |2l 

In this section we study the dipolar motion of ideal bosonic 
and fermionic gases trapped in the combined lattice and har- 
monic potentials. We start by writing an expression for the 
evolution of the center of mass of an ideal gas with general 
quantum statistics, and then we use this expression to study 
the dipole oscillations for bosonic and fermionic systems. The 
simplicity of the noninteracting treatment allows us to derive 
analytic equations for the dipole dynamics for both statistics. 

Later on, in section lTVl we show how the knowledge of the 
bosonic and fermionic ideal gas dynamics can be useful in 
describing the dynamics of the interacting bosonic system for 
a large range of parameters of the trapping potentials. 
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A. Ideal gas dynamics 

Consider an ideal gas of N atoms loaded in the ground state 
of an optical lattice plus a quadratic potential initially dis- 
placed from the trap center by 5 lattice sites. The initial state 
of the gas is described by a mixed ensemble state with mean 
occupation numbers determined by the appropriate quantum 
statistics 



n 

™n = fi„-l , ( 27 ) 



If the initial displacement of the atomic cloud is small, 
25 <C n c , and the lattice is not very deep (q ^> 1), local- 
ized eigenstates are initially not populated. Then, only low- 
energy states are relevant for the dynamics and the latter can 
be modeled by utilizing the asymptotic expansions derived in 
Sec. Ill Bl To simplify the calculations, we use the harmonic 
oscillator approximation for the eigenmodes (Ea.( ll6» . and in- 
clude up to the quadratic corrections in n in the eigenenergies, 
which corresponds to keep the first three terms of Ea. (ll4» . 
Even though this treatment is not exact, we found that it prop- 
erly accounts for the period and amplitudes of the center of 
mass oscillations for small trap displacement. After some al- 
gebra it is possible to show that the time evolution of the center 
of mass is given by 



where T is the initial temperature of the system, fc^ is the 
Boltzmann constant, /i is the chemical potential which fixes 
the total number of particles to iV, ^ rT„ = N, and the posi- 
tive and negative signs in Eq.( I27> are for fermions or bosons, 
respectively. 

The time evolution of the center of mass of the gas is dic- 
tated by the ensemble average 



(28) 



with 



( Xn (t)) = aJ^l^^e-^-^Eif^f: 

k,l 



ifc) f (l)* 



(n) f (k) 
j ' 



(29) 



where the quantities / . and the energies correspond to 
eigenvalues and eigenenergies of the undisplaced system. The 

(n) 

coefficients c k are given by the projection of the n excited 
displaced eigenstate onto the k excited undisplaced one. 

Once the E n and /j™' are known, the center of mass evo- 
lution can be calculated. In the following we discuss the zero 
temperature dynamics for the ideal bosonic and fermionic sys- 
tems. 



1. Bosonic system 

At zero temperature the bosons are Bose condensed and 
n n = NS n Q, where 5 n o is the Rronecker delta function. The 
center of mass motion is then given by 



it)) = a£ «f c 



( )J°) o -i(E h -E t )t/H 



(As) f (l)» 
3 



k,l 



JO) 



(o) f (k) 
j ' 



(30) 



= aSe 



(sn) 



cos 0J„t — 



'OA 
■ sm l — I 



(31) 



with fad* = Hlo* - 0/4. 

In Fig. 0]we plot the average center of mass position in 
units of the lattice spacing a as a function of time for an ideal 
bosonic system of atoms with 87 Rb mass. The solid line is ob- 
tained by numerically solving the tight binding Schrodinger 
equation, Eq.Q, while the dotted line is the analytical so- 
lution Eq.(l3U. For the plot we used V = 7 AEr, cut — 
2ir x 60Hz and 5 = 3. The time is shown in units of 
T Q = 2ir/ui*, a characteristic time scale. The two solutions 
exhibit very good agreement for the times shown. 

The modulation of the dipole oscillations predicted by 
Ea. (l31> can be observed in the plot. At early times, t <IC H/tl, 
the amplitude decreases exponentially as exp(— T t 2 ), with 

r o = (^ s 2^ h ) ' an d tne frequency is shifted from u* by 

|| (1 — A-), The initial decay does not correspond to real 
damping in a dissipative sense, as in a closed system the en- 
ergy is conserved. The decay is just an initial modulation and 
after some time revivals must be observed. Because in the 
large q limit ui'q <C fl/h, the revival time is approximately 
given by ih/fl. 

It is a general result that the dipole oscillations of a harmon- 
ically confined gas in absence of the lattice are undamped. 
The undamped behavior holds independently of the temper- 
ature, quantum statistics and interaction effects (generalized 
Kohn theorem 1 26]). Equation ( 13 11 shows how this result does 
not apply when the optical lattice is present even for an ideal 
Bose gas. 

Recent experimental developments have opened the possi- 
bility to create a non-interacting gas for any given strength of 
the trapping potentials l27l l28ll . The techniques utilize Fes- 
hbach resonances for tuning the atomic scattering length to 
zero. These developments should allow for the experimental 
observation of the modulation of the dipole oscillations of an 
ideal gas predicted in this section. 
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2. Fermionic system 
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FIG. 4: Center of mass motion in lattice units as a function of time 
for an ideal bosonic gas. In the plot, V = IAEr, lot = 2n x QQHz 
and 5 = 3. The time has been rescaled by T = 2n /lu* . The solid 
and dotted lines are the numerical and analytical solution Eq. J3 II . 
respectively. 
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FIG. 5: Center of mass motion in lattice units as a function of time 
for N = 15 fermions. Here lot = 2i\ x 2QHz and V a = 7AEr and 
5 = 3. The time is in units of T = 2tt/ui* . The solid and dotted 
lines are the numerical and analytical solution Eq. 1341 . respectively. 



J 



At zero temperature the Pauli exclusion principle forces 
fermions to occupy the lowest N eigenmodes, and therefore 
n n = 1 for < n < N — 1 and zero elsewhere. The oc- 
cupation of the first N displaced modes makes the condition 
of occupying only low-energy eigenstates of the undisplaced 
potentials more restrictive than in the bosonic case. Neverthe- 
less, if the initial displacement, lattice depth and atom number 
are chosen such that only LE eigenstates are initially popu- 



lated, 



„(JV-1)|2 



<C 1, it is possible to derive simple 



analytic expressions for the dipole dynamics. This is the fo- 
cus of the remainder of this section. Population of localized 
states considerably complicates the system's dynamics and a 
numerical analysis is therefore required. This is postponed to 
Sec lIVB2l 



When only LE undisplaced eigenstates are occupied, as 
explained for the bosonic system, to a good approximation 
the eigenmodes can be assumed to be the harmonic oscilla- 
tor eigenstates and only corrections quadratic in the quantum 
number n are relevant in Ea. (ll4> . After some algebra, the 
above approximations yield the following expression for the 
time evolution of the center of mass: 



JV-l 



n=0 



(x n (t)) = aSKe exp 



x n (t) 




( x (t) - i) 2k x(f) n - 1 - k ((n + l)x(t) + k - n) 



(2fc)!!(fc + 1)! 
where X (t) = exp(- 40t /( 4fi » and fiw D * = hw* - Vt/A. 



]](n + l-s) 



(£) (i-xW 



s=l 



(2n)\\ 



(32) 
(33) 

n 

-(34) 



r 



The parameter \(t) takes into account the quadratic cor- tice. In the limit \(t) — > 1, x n (t) — * 1 Vrt and therefore the 
rections to the harmonic oscillator energies. The corrections amplitude of the dipole oscillations remains constant in time, 
are proportional to f2/4, and due to the presence of the lat- (x) — dS cos(uj*t), as predicted by Kohn theorem. The cor- 
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FIG. 6: Center of mass motion in lattice units as a function of time for some initially occupied modes. The modes are labeled by the quantum 
number n. As in Fig. [5] the parameters are lot = 2-k x 20Hz, V = 7AEr, N = 15, 6 = 3 and the time is in units of T a — 2-k/ui* . The 
solid line corresponds to the numerical solution and the dotted line to the solution given by Eq. 1341 . When the center of mass evolution of all 
the different modes is added, from n = to n = N — 1 one recovers the total center of mass evolution shown in Fig. [3] 



rections quadratic in n cause the modulation of center of mass 
oscillations. 

The modulation is caused not only by the over- 
all envelope generated by the exponential term 
cxp (— S 2 (l — x(t))/(2a| )), which was also present in 
the bosonic case, but mainly from the interference created 
by the different evolution of the N average positions (x n ) in 
the sum Ea.( l32t . The latter induces a fast initial decay of the 
amplitude of the dipole oscillations. 

In Fig|5]we plot the center of mass motion of the fermionic 
gas composed of N — 15 atoms with the mass of 87 Rb. The 
solid and dotted lines correspond to the numerical and analytic 



solutions, respectively. Here the depth of the optical lattice is 
7.4Er, and lot = 2tt x 20Hz. The amplitude of oscillation 
shows a rapid decay in time. The analytic solution captures 
the overall qualitative behavior of the numerical curve. Nev- 
ertheless, only at short times the agreement is quantitatively 
good. Population of eigenstates which are not fully harmonic 
in character is responsible for the disagreement at later times. 
This effect is particularly relevant for the evolution of the dis- 
placed states with larger quantum number, as explicitly shown 
in Fig. [S]where the time evolution of some displaced modes is 
plotted. Again, the solid line is the numerical solution and the 
dotted line is the analytic one. For the lowest energy modes, 
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n = and n = 3, the agreement between the two curves is 
almost perfect. For the higher energy modes n = 6 and n = 9 
the analytic solution is underdamped and overestimates the 
collapse time. 

Interestingly, the dynamics of the displaced excited modes 
exhibits an initial growth of the amplitude. This behavior 
is a pure quantum mechanical phenomenon due to the 
constructive interference between the different phases of the 
undisplaced eigenmodes during the evolution. We explicitly 
checked for energy conservation during the time evolution. 
The amplitude increase is captured by the analytic solution 
and it allowed us to show that the growth happens only 
when the ratio between the initial displacement 5 and a^o 
is less than one. While such a behavior is not observable in 
the evolution of a fermionic cloud, as the observable is the 
center of mass position summed over all initially populated 
modes (x(t)), the experimental observation of growth for an 
individual mode may be possible if an ideal bosonic gas is 
initially loaded in a particular excited state, and then suddenly 
displaced. 

As described above, the evolution of LE modes can be 
handled analytically. On the other hand, when high-energy 
eigenmodes are populated the dynamics is much more com- 
plicated. Nevertheless, there is another simple limiting case 
that can be actually solved. This corresponds to the case 
when the displacement is large enough or the lattice deep 
enough that the displaced cloud has non-vanishing projec- 
tion amplitudes only onto high-energy undisplaced modes 
which can be roughly approximated by position eigenstates, 
j r, r— ^ (5j. r ±5-j. r )/y/2. Then, one finds that (x n ) sa aS 
for all n and thus (x) ~ a5. That is, when only high-energy 
eigenmodes are populated the dynamics is completely over- 
damped and the cloud tends to remain frozen at the initial 
displaced position. We show later on in Sec. IIV B 21 where 
we treat interacting atoms, that in the so-called Mott insulator 
regime most populated modes are actually localized, and this 
kind of overdamped behavior characterizes the dipole dynam- 
ics. 



IV. MANY-BODY SYSTEM 

A. Spectrum of the BH-Hamiltonian in presence of an external 
quadratic potential 

The Bose-Hubbard (BH) Hamiltonian describes the sys- 
tem's dynamics when the lattice is loaded such that only the 
lowest vibrational level of each lattice site is occupied 1 29] 

U 



H 



BH 



E 



sir 



J{a]a 



a j+i a . 



) + ^{fij - 1) 
(35) 

Here dj(aA is the bosonic annihilation(creation) operator of 



a particle at site j and hj = 



SI and J are defined as 



in Eqs. @ and (jSJi, while U is an on-site interaction energy 



given by U = 



J dx\w(xo)\ 4 , with a s the s-wave scat- 



at the origin. The quantity U is the energy cost for having two 
atoms at the same lattice site. The tunneling rate decreases for 
sinusoidal lattices with the axial lattice depth V Q as 



J = A 



E R 



exp 




(36) 



where the numerically obtained constants are A = 1.397, 
B = 1.051 and C = 2.121. The interaction energy increases 
with V as 



Er 



1/4 



(37) 



where is a dimensionless constant proportional to a s . 
In current experiments, the one-dimensional lattice is ob- 
tained by tightly confining in two directions atoms loaded 
in a three-dimensional lattice. In this case (3 = 
4%/27r(a s /A)(Vj_/-E_R) 1/2 , where V± is the depth of the lattice 
in the transverse directions 1 30, 31]. The parameter 7 = U/J 
therefore increases as a function of the axial lattice depth as 



1{V ) 



I 
A 



Vo_ 

Er 



1/4-fl 



exp 




(38) 



tering length and wq [x) the first-band Wannier state centered 



In the absence of the external quadratic potential, the 
bosonic spectrum is fully characterized by the ratio between 
the interaction and kinetic energies 7 and the filling factor 
N/M, where M is the number of lattice sites |32]. For 7 <C 1 
and any N/M ratio the system is weakly interacting and su- 
perfluid. For 7 » 1 and M > N the system fermionizes 
to minimize the inter-particle repulsion. In this regime the 
bosonic energy spectrum mimics the fermionic one, the corre- 
spondence being exact in the limit of infinitely strong interac- 
tions. In particular, in a lattice model the onset of fermioniza- 
tion is characterized by suppression of multiple particle occu- 
pancy of single sites. This implies that fermionization occurs 
for eigenstates whose energy is lower than the interaction en- 
ergy U. If M > N there are M\/{N\{M — N)\) fermionized 
eigenmodes and the dynamics at energies much lower than U 
can be accounted for by using these states only. On the other 
hand, if the lattice is commensurately filled, N — M, there 
is only one fermionized eigenstate, and it corresponds to the 
ground state. All excited states have at least one multiply oc- 
cupied site, and therefore excitations are not fermionized. The 
ground state corresponds to the Mott state with a single par- 
ticle per site and reduced number fluctuations. The transition 
from the superfluid to the Mott state is a quantum phase tran- 
sition, and the critical point for one-dimensional unit filled 
lattices is j c ~ 4.65 132113 311. 

In the presence of the quadratic trap the spectrum is deter- 
mined by an interplay of U, J, SI and N. In trapped systems 
the notion of lattice commensurability becomes meaningless 
because the size of the wave-function is explicitly determined 
by these parameters. As a consequence, for any value of 
the ground state can be made to be a Mott insulator with one 
atom per site at the trap center by an appropriate choice of 
U, J and Si [34], and the lowest energy modes can always be 
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FIG. 7: Energy spectra as a function of the depth of the axial optical 
lattice Vo. The continuous, dotted and dashed lines correspond to 
TV = 5 interacting bosons, non-interacting bosons and fermions, re- 
spectively. The dashed-dotted line corresponds to hcu* . The horizon- 
tal axis on the top of the figure is 7 = U / J, and is only meaningful 
for interacting bosons. For each energy spectrum, the corresponding 
ground-state energy Eo has been subtracted. 



made to be fermionized in the large U limit. The purpose of 
this section is to characterize fermionization and localization 
of the many-body wave-function when both the quadratic and 
periodic potentials are present, by relating the occurrence 
of the different regimes to changes in the spectrum at low 
energies. 

We performed exact diagonalizations of the BH- 
Hamiltonian for N = 5 particles and M — 19 sites in 
presence of a quadratic trap of frequency lot — 2tt x 150 Hz. 
For the calculations, we chose 87 Rb atoms with scattering 
length a s = 5.31 nm, a lattice constant a = 405 nm, and 
therefore ~ 0.0046-Er. We fixed the transverse lattice 
confinement to V± ~ 25.5 E r and varied the depth of the 
optical lattice V in the parallel direction from 2 Er to about 
17 En- For these lattice depths, the energies U and J both 
vary so that their ratio 7 increases from 3.3 to about 150. 
Due to the changes in J, the ratio q = characterizing 
the single-particle solutions decreases from 130 to 4 with in- 
creasing V . The effective harmonic-oscillator energy spacing 
tuo* decreases approximately from 0.0525 Er to 0.0092 Er. 
We used these parameters because they are experimentally 
feasible and fulfill the condition U - fl((N - l)/2) 2 > 
for the entire range of the trapping potentials. Later on we 
discuss that, for deep enough lattices, fulfillment of the last 
inequality ensures the existence of an energy range in which 
eigenmodes are fermionized. 

Figure Q shows the lowest eigenergies of the BH- 
Hamiltonian as a function of V . The continuous line is the 
exact solution for N = 5 interacting bosons. The dotted 
and dashed lines are the exact spectra for 5 non-interacting 
bosons and fermions, respectively. They have the same mass 
and are trapped in the same potentials as the interacting 



bosons. Their spectra are shown for comparison purposes. 
For each spectrum the energy Eo of the ground state has been 
subtracted. 

In the absence of the optical lattice, V = 0, the energy 
difference AEq = E\ — Eo, or energy spacing, between the 
first excited and ground state equals the harmonic oscillator 
level spacing Itiot, independent of statistics and interaction 
strength. Figure show that this is no longer the case in the 
presence of the lattice. 

The dependence of AEo on statistics is evident in the 
plot, as the level spacing is different for ideal bosons and 
fermions. In particula r, the energy spacing for bosons is only 
shifted from fko* = \/4£!J (dashed-dotted line) by an amount 
which is almost constant for all lattice depths, while AEo f° r 
fermions clearly deviates from Hui*, especially for deep lat- 
tices. The behavior of AEo m the two cases can be understood 
by using the asymptotic solutions of the single-particle prob- 
lem. For ideal bosons AEo is equal to the energy difference 
between the first-excited and ground single-particle eigenen- 
ergies. The ratio q used for the plots is such that the critical 
value n c of Eq. j23t is always larger than 2, and therefore the 
ground and first-excited eigenenergies are well described by 
Eq. dl4> . The calculation of the energy difference using this 
equation yields AEo ~ ^* — 0/4. For fermions, AEo is 
equal to the difference between the energies of the n = N and 
n = N — 1 single-particle excited states. For V < 9.6Er, the 
critical value n c is larger than TV, and therefore the energies 
of the n — N and n = N — 1 single-particle excited states 
are also well described by Eq. (114-1 . Then, AEo f° r fermions 
is smaller than for bosons because lattice corrections are more 
important for higher quantum numbers, and have all negative 
sign. On the other hand, for V > 9.6Er, n c is smaller than 
N — 1 and the energies of the n = N and n = N — 1 single- 
particle excited states are described by Eq. J20I . The transition 
of the single-particle eigenmodes at the Fermi level from LE 
to HE around Vq = 9.6Er is signaled by the minimum of 
AEo for fermions. In general, an estimate for the value of J 
at which the minimum takes place is 

J ~ n((7V — l)/2) 2 /2. (39) 

This value is obtained by equating the Fermi energy En-i, 
which is of order Q((N - 1)/2) 2 from Eq. @D), to E 7lc , which 
is approximately 2J. The transition of the single-particle 
eigenmode at the Fermi level from LE to HE is also con- 
nected to the formation of a region of particle localization 
at the trap center in the many-body density profile. As ex- 
plained in 1 13], when £jv-i is equal to E Ua the on-site den- 
sity in the central site of the trap approaches 1 with reduced 
fluctuations. For Vq > 9.6Er, Fig. shows that AEo ap- 
proaches an asymptotic value flN, value that can be derived 
from Eq. ( I20> . When the asymptotic value £IN is reached, 
most single-particle states below the Fermi level are localized, 
and this yields a many-body density profile with N unit-filled 
lattice sites at the trap center. 

The dependence of the first excitation energy on interac- 
tions can also be seen in Fig.0 In fact, by comparing AEo 
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for the interacting bosons to the value of AEq for ideal bosons 
and fermions, three different regimes can be considered: 1 < 
7 < 10, 10 < 7 < 30 and 7 > 30. These regimes cor- 
respond to the intermediate, fermionized-non-localized and 
fermionized-Mott regimes, respectively. The weakly inter- 
acting regime 7 < 1 is only reached for V -C 2Er for 
our choice of atoms and trapping potentials. For such lat- 
tice depths the tight-binding approximation is not valid, and 
therefore we do not show the spectra for this regime in Fig. [7] 
In the following we discuss the main features of the different 
regimes focusing on the connection to the ideal bosonic and 
fermionic systems. 

• For 7 < 1 the interacting bosonic system is in the weakly 
interacting regime. In this regime the first excitation energy 
is almost the same as the ideal bosonic one. Most atoms 
are Bose-condensed, interaction-induced correlations can be 
treated as a small perturbation, and the spectrum can be shown 
to be well reproduced by utilizing Bogoliubov theory [ 35 ] . 

• For 1 < 7 < 10, the system is in the intermediate regime, 
where Ai^o for interacting bosons deviates from the ideal 
bosonic energy spacing and approaches the ideal fermionic 
one. Indeed, FigQ shows that for Vo < 4Er, AEq for the 
interacting bosons lies closer to the ideal bosonic energy spac- 
ing, while for Vq > 4 it lies closer to the ideal fermionic one. 
In the presence of the optical lattice 7 increases exponentially 
with V ot and therefore the intermediate regime occurs for a 
relatively small range of accessible trapping potentials, here 
for 1 < V /E R < 5.5. 

• For 7 > 10 the interacting spectrum approaches the ideal 
Fermi spectrum and the system is in the fermionized regime. 
The numerical solutions show that the energy difference be- 
tween the energy spectra of interacting bosons and fermions 
is of the order of J 2 /U and slightly increases for larger fre- 
quencies of the quadratic trap. 

In general, fermionization in the presence of the external 
quadratic potential occurs for N < M when the two following 
inequalities are satisfied 



7 = U/ J > 1, U > n((N - l)/2) 2 



(40) 



While the first inequality is the same as for homogeneous lat- 
tices and relates to the building of particle correlations, the 
second inequality is specific to the trapped case and relates to 
suppression of double particle occupancy of single sites. If the 
interaction energy U is larger than the largest trapping energy, 
which corresponds to trapping an atom at position (N — l)/2, 
it is energetically favorable to have at most one atom per well. 
The average on-site occupation is therefore less or equal to 
one. 

The second inequality in Eq. (I40i poses some limitations 
on the choice of possible 51 and U for a given number of 
trapped particles. For our choice of the trapping potentials, 
this inequality is satisfied for any V a . Indeed, this is not 
an unrealistic assumption. In recent experiments with 87 Rb 
atoms, an array of fermionized gases has been created with 
at most 18 atoms per tube 0]. For such N, the condition 
U > fl{{N - l)/2) 2 can be fulfilled for many different 
choices of experimentally feasible trapping potentials. 
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FIG. 8: Density profiles rij = (hj) as a function of the lattice site j 
for different lattice depths. The dash-dotted, dashed and dotted lines 
correspond to interacting bosons, while the boxes, dots and trian- 
gles correspond to ideal fermions for V /Er = 7, 12 and 15 lattice 
depths, respectively. 
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FIG. 9: Number fluctuations An-i 



i] = (n?) — {fij) 2 as a function 
of the lattice site j for different lattice depths. Conventions and pa- 
rameters are the same as in Fig|8| 



• For 7 > 30 the system enters the fermionized-Moff 
regime. In fact, Fig.0shows that for 7 > 30 the energy spac- 
ing for the interacting bosons (and also for ideal fermions) 
begins to increase until it reaches an asymptotic value at 
7 « 150. 

The approaching of the asymptotic value signals the forma- 
tion of a localized many-body state for the interacting bosons, 
the so-called Mott insulator state. The relevant relation be- 
tween N, fl and J for the formation of an extended core of 
unit-filled sites at the trap center with fluctuations mainly at 
the outermost occupied sites is 



flN > J. 



(41) 



This is explained as £IN is the energy cost for moving a par- 
ticle from position \ (N - 1)/2| to position \ (N + 1)/2| at the 
borders of the occupied lattice and it is the lowest excitation 
energy deep in the Mott state. 

Figure shows that for 7 « 150 the energies of the first 
four excited states become degenerate. This degeneracy 
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occurs because deep in the Mott regime the energy required 
to shift all the atoms of one lattice site to the right or left 
is the same as the energy required for moving an atom 
from site ±(N - l)/2 to site ±{N + l)/2. For 7 ps 150, 
J is approximately fl, and therefore the tunneling energy 
is barely sufficient to overcome the potential energy cost 
D, for moving a particle from the central site of the trap 
to one of its neighbors. In the single-particle picture, when 
J ps fl all the single-particle states below Ejf-i are localized. 

In order to better understand the formation of the Mott in- 
sulator, in Figs. [3]and[9]the on-site particle number Uj = (fij) 

and fluctuations Arij = \J (n?) — {hj) 2 are plotted as a func- 
tion of the lattice site j, for different lattice depths V a . The 
lines and symbols are the results for interacting bosons and 
ideal fermions, respectively. All the V Q -values are such that 
the interacting bosonic system is fermionized. This is mir- 
rored by the overall good agreement between lines and sym- 
bols for all the curves. The plots show that for V = 7Er{^ = 
15, dashed-dotted line), the largest average particle occupa- 
tion is no fa 0.7, and number fluctuations are of the order 
of 0.5 in the central 9 sites, while for V — 12Er(^ = 54, 
dashed line), n approaches 1 in the central 3 sites, and fluc- 
tuations at the trap center drop to a value Arij pa 0.1. The 
sharp drop in particle fluctuations clearly signals the localiza- 
tion at the trap center, and is in agreement with J < UN for 
7 = 54. For V = 15E R ('y = 104, dotted line), the Mott state 
is formed, as the mean particle number in the five central sites 
is one, with nearly no fluctuations. Fluctuations are larger at 
lattice sites far from the center, and due to the tunneling of 
particles to unoccupied sites. 

Because of the small number of atoms that we use in the cal- 
culations, Q((N — l)/2) 2 /2 and QN are of the same order of 
magnitude. It is therefore not possible to clearly distinguish 
the value of J for which the on-site density at the trap center 
approaches one from the value of J for which the Mott state 
is fully formed at the trap center, with reduced number fluc- 
tuations in the central N sites. Preliminary results obtained 
with a Quantum Monte-Carlo code, Worm Algorithm |41], 
confirm the existence of these two distinct parameter regions 
when more atoms are considered, and therefore the usefulness 
of both the energy scales Q((N - l)/2) 2 /2 and QN for inter- 
acting bosons. These results will be published elsewhere. 



B. Many-body dynamics 

In this section we study the temporal dipole dynamics of an 
interacting bosonic system composed of 5 atoms trapped in 
a combined quadratic plus periodic confinement. The role of 
interactions on the effective damping of the dipole oscillations 
is studied by means of exact diagonalization of the Hamilto- 
nian. As discussed when dealing with the ideal gas dynamics, 
such damping is effective because it is due to dephasing and 
does not have a dissipative character. 

Assuming the system is initially at T — 0, the evolution of 
the center of mass is given by: 



(x(t)) = Y, A ike luJlkt (42) 

l,k 

A lk = aY,3{<t>iM<t>k)CtC k , (43) 
3 

with hid it = Ei — Ek and where Ei and \<f>i) are eigenvalues 
and eigenmodes of the Z?//-Hamiltonian. The coefficients C; 
are the projections of the initial displaced ground state |<^(0)) 
onto the eigenfunctions {| cf>)i} of the undisplaced Hamilto- 
nian,G = (0/1^(0)). 

For the exact evolution, the ground-state |y(0)) is calcu- 
lated by shifting the center of the quadratic trap by S lattice 
sites. The number of wells M is 19 for all simulations, which 
fixes the size of the Hilbert space to 33649. In the time prop- 
agation we only keep those eigenstates whose coefficients C; 
are such that |C/| 2 > 10~ 3 . The typical number of states 
that fulfil this requirement is about 100. The accuracy of the 
truncation of the Hilbert space during the time propagation 
has been checked by increasing the number of retained states, 
finding no appreciable changes in the dynamics. 

We are interested in the dynamics both when a Mott insu- 
lator state is not and is present at the trap center. These two 
cases are discussed in Secs lIV BTl and lTV B 21 respectively. In 
particular, in Sec. II V B li the interaction strength U is varied, 
while the ratio q, specifying the ideal gas dynamics, is large 
and constant. For the chosen values of q and N, J > UN 
and for increasing U the system fermionizes without forming 
a Mott insulator at the trap center. In Sec. IIV B 21 the dynam- 
ics of systems that do exhibit a Mott insulator in the large U 
limit is analyzed. In this case, J and U are simultaneously 
varied by increasing the axial optical lattice depth. 

1. Non-localized dynamics 

In the absence of the optical lattice, the equations of mo- 
tion for the center of mass are decoupled from those of the 
relative coordinates. As only the latter are affected by inter- 
actions, all modes excited during the collective oscillations 
have the harmonic oscillator energy spacing fiuy, and there- 
fore (x(t)} = (x(0)} cos(w T £). 

When the lattice is present, the equations of motion for 
the center of mass and relative coordinates are coupled, and 
thus the many-body dynamics is interaction dependent. In 
this section we fix lot — 2tt x 100Hz and V = 7Er, and 
study the role of interactions in the many-atom dynamics by 
varying 7 from to 200, for constant q = 77. This can be 
experimentally realized by tuning the scattering length of the 
system by means of Feshbach resonances. In the following, 
we analyze the weakly interacting, intermediate and strongly 
interacting regimes separately. 

Weakly interacting regime: 7 < 1 

In order to study the role of interactions in the weakly in- 
teracting regime, in Fig. ^| the effective damping constant of 
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dipole oscillations T is shown as a function of 7. The damping 
constant was calculated by fitting the first 10 center of mass 
oscillations to the ansatz (x(t)) = aScxp (— Ft 2 ) cos(u>t), 
where V and uj are fitting parameters. This ansatz is cho- 
sen in analogy to the non-interacting model. The effective 
damping T is in units of T = 8 2 £l 2 / (8haho) 2 , which is the 
damping constant predicted by Eq. The solid and dotted 
lines correspond to T as calculated by means of exact diag- 
onalizations and by numerically evolving the following Dis- 
crete Non-Linear Schrddinger Equation (DNLSE) for the am- 
plitudes {Zj} 

dz- 

ih-gj- = -J{zj+i + Zj-i) + VLj 2 Zj + U\ Zj \ 2 z j: (44) 

respectively. Ea.J44> has been obtained by replacing the field 
operator ctj with the c-number Zj (t) in the Heisenberg equa- 
tion of motion for a,j. Such replacement is justified for 7 <§; 1 
as the many-body state is almost a product over identical 
single-particle wave-functions. The amplitudes {zj} satisfy 
the normalization condition J^j \ z j\ 2 = N. The initial state 
used in the evolution of the DNLSE was found by numerically 
solving for the ground state of Eq. (1441 . displaced by 8 = 2 
lattice sites. 

1 pi 1 1 1 1 n 
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FIG. 10: Effective damping constant of the dipole oscillations as a 
function of 7 for a system in the weakly interacting regime. Here, 
q w 77, lot = 100 Hz and r o = (SQ/Sha ho ) 2 . 

In Fig. El| the continuous and dotted lines overlap for 
7 < 0.05, and show a decrease in the damping constant with 
increasing interaction strength in the whole range 7 < 0.2. 
For values of 7 > 0.05 the mean field and exact solutions start 
to disagree. While the exact solution has a minimum around 
7 « 0.2, and then grows for larger 7 values, the mean field 
curve decreases monotonically to zero. 

The fact that the mean-field solution decreases to zero for 
increasing interactions is explained by noting that when inter- 
action effects become important the density profile acquires 
the form of an inverted parabola, or Thomas-Fermi profile, 
\zJ F {t)\ 2 « n {(j - (x(0))/af - R 2 TF ) /U, (x(0)> = aS. 
Substitution of the Thomas-Fermi profile in Eq. (1441 leads to 
the exact cancellation of the quadratic potential, and thus, in 
the frame co-moving with the atomic cloud, the atoms feel an 
effective linear potential. The spectrum of a linear plus peri- 
odic potential is known to be equally spaced |36], and there- 
fore no damping due to dephasing is expected. 



It is important to note that the mean-field undamped oscil- 
lation occurs only in a parameter regime far from dynamical 
instabilities. In fact, as shown in previous theoretical and 
experimental studies when the initial displace- 

ment is large enough to populate states above half of the 
lattice band-width, mean-field dynamical instabilities induce 
a chaotic dipole dynamics. In the framework of this paper, 
this critical displacement corresponds to 5 w n c /2. For 

6 > n c /2, the initial ground-state has a significant overlap 
with localized eigenmodes of the undisplaced system, which 
are therefore populated during the dipolar dynamics, causing 
damping. The importance of the population of these modes 
is enhanced by the non-linear term, which causes an abrupt 
suppression of the center of mass oscillations at the critical 
point in the mean-field solution. 

While the mean-field anaysis accounts for the decrease in 
the damping constant, Fig.^|shows that it is not accurate for 

7 > 0.2. This is due to the fact that the mean-field analysis 
completely neglects interaction-induced correlations. These 
correlations are responsible for the quantum depletion of the 
condensate, which causes some atoms to be excited to higher- 
energy single-particle eigenmodes which are more affected by 
the lattice. To show this effect in a quantitative manner, in 
Fig[TJwe plot the probability distribution of the frequencies 
wz/c, Eq. d42i . for some values of 7. In the histograms, the 
height of a bar-chart centered at a given frequency u is the 
occurrence probability of uiik- The probability is proportional 
to the normalized sum over all the weight factors \Aih | whose 
frequency lies between toik — 0.0025 and u>ik + 0.0025. The 
frequencies are in units of the effective harmonic oscillator 
frequency to* of Eq. dl7l . This approach is similar to the one 
used in Ref. |42], where the strength function is used to study 
the collective dynamics induced by mono- and dipolar excita- 
tions. 

The histogram for the case 7 = shows a frequency distri- 
bution with most frequencies in the interval 0.85a;* < u>ik < 
0.96a;*. In particular, two large peaks are observed in the 
range 0.9a;* < u>ik < 0.96a;*. This is to be compared to 
the case where the lattice is not present, and a single peak at 
to* is expected. The observed frequency spread is due to the 
modification of the harmonic oscillator spectrum introduced 
by the lattice and is responsible for the observed damping in 
the ideal bosonic gas, as explained in detail in Sec. MI ATI 

Figure \H\ also shows that for 7 = 0.3 the system has a 
narrower frequency distribution. In this case approximately 
100% of the frequencies lie in the interval 0.9a;* < ujik < 
0.96a;*. The frequency narrowing from 7 = to 7 = 0.3 is 
consistent with the decrease in the damping constant shown 
in Fig. [TO] For larger values of 7, 7 = 1 and 2, some modes 
with frequency smaller than 0.7a;* and larger than l.lui* start 
to contribute to the collective dynamics. Population of such 
modes is related to the depletion of the condensate and sig- 
nal the increased importance of quantum fluctuations in the 
system. 

Finally, we note that for our choice of il, N and J, 
f2((JV — l)/2) 2 /J is approximately 0.2, which is the value 
of 7 = U/J at which the mean-field and exact solutions 
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FIG. 11: Probability distribution of frequencies for different values of 7, with q « 77 and ut = 100 Hz. The frequencies 0% are in units of 



begin to disagree. This suggests that the fulfillment of the 
second inequality in Eq. EJ, U > Q((N ~ l)/2) 2 , is re- 
lated to the failure of the mean-field approach, even for 7 < 1 . 

Intermediate Regime: 1 < 7 < 10 

In Fig. the numerically obtained damping constant is 
shown for 7 > 1. In this parameter regime we find that the 
function ad cxp(— Tt 2 ) cos(wi) does not provide a good fit to 
the center of mass evolution. Instead, we find that a better fit- 
ting ansatz is given by a<5 exp(— Tt) cos(wi) . In the plot, the 
damping constant is normalized to T\, which is the damping 
rate at 7 = 1. 

We observe that for 2 < 7 < 5 the damping is almost 
constant. This is consistent with the fact that by inspection 
the spectrum of excited frequencies has a similar shape and 



width in the entire transition region. An example of such a 
frequency distribution is given in Fig. ^2 f° r 7 = 4. The 
dominant peak is around uiik ~ 0.9w*, while multiple peaks 
are noticeable between 0.7oj* < oJik < 0.92a;*. The overall 
envelope of the distribution has a long tail, as opposed to the 
case 7=1 where all the weight is roughly concentrated in 
just two frequencies. 

The increased importance of the tails of the distribution 
for 7 > 1 qualitatively accounts for the transition from an 
exponential decay quadratic in time towards an exponential 
decay which is linear in time. In fact, for 7 < 1 the proba- 
bility distribution of frequencies may be fitted by a Gaussian, 
while for 7 > 1 a better fit is provided by a Lorentzian-like 
profile. The Fourier transforms of such distributions give 
precisely the observed functional form of the decay of the 
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FIG. 12: Effective damping rate of the dipole oscillations T as a 
function of 7 in the intermediate and strongly correlated regimes. 
Here, q w 77 and lot = 100 Hz. The damping rate has been rescaled 
such that r(7 = 1) = 1. The large and small dots are for interacting 
bosons and non-interacting fermions, respectively. The continuous 
line is a guide for the eye. The dashed line is the best-fit curve Y = 
10.13 exp(— 16. 47 -1 ' 25 ) to the exact damping rate for interacting 
bosons. 



dipole oscillations. 

Strongly interacting regime 7 > 10 

In Fig. E|f° r 7 > 10, the damping rate is shown to 
rapidly increase and approach a finite asymptotic value which 
is depicted in the plot by a dotted line. This asymptotic 
value Too corresponds to the damping rate calculated for an 
ideal fermionic gas. The fermionic damping rate is con- 
stant because here J and £1 are kept constant while U in- 
creases. The tendency to approach the fermionic damping 
rate as 7 increases is a consequence of fermionization of the 
bosonic wave-functions for 7 > 1, Eq. ( 140b . Numerically 
we find that the damping rate approaches exponentially, 
r(7) = Too exp(— aj a ), with a best-fit exponent a of order 
— 1. The fitting curve is shown in the plot with a dashed line. 

The dipole dynamics of the bosonic and fermionic systems 
are explicitly compared in Fig.^J where we plot the first 10 
oscillations of the center of mass after the sudden displace- 
ment of the trap. In the figure, the dashed line and the dots 
are the bosonic dynamics for 7 = 13.6 and 100, the solid 
line is the fermionic evolution, and the crosses are the analyti- 
cal approximation to the fermionic evolution Eq. ( I34> . respec- 
tively. Consistent with Fig.^J we observe that for increasing 
7 the decrease of the amplitude of oscillation for the bosons 
resembles more and more the one for fermions. In particu- 
lar, for 7 = 100, the curves for the interacting bosons and 
ideal fermions nearly overlap. The distribution of excited fre- 
quencies for 7 = 100 is shown in Fig. ^2 The frequency 
distribution is broad and centered around u> w 0.75oA Also 
in the inset small peaks are shown to appear in the range 
l.6u* < u> < 3oj* (notice the different scale in the inset). The 
broad distribution is due to the population of single-particle 
states that are not harmonic in character. For the value of q 
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FIG. 13: Center of mass position in lattice units as a function of time 
for interacting bosons, for 7 = 13.6(dashed line) and 7 = 100(dots). 
The solid line is the center of mass position for ideal fermions, while 
the crosses are the analytical approximation to the fermionic solution 
Eq.j34}. Here, q « 77 and co T = 2tt x 100 Hz and T a = 2it/lo*. 



used for the calculations no single-particle localized modes 
are occupied in the ground state before the trap displacement. 
After the displacement about 90 percent of the atoms occupy 
non localized single-particle modes. The phase mixing be- 
tween these modes accounts for most of the observed damp- 
ing. The remaining 10 percent occupy localized states and the 
population of these modes is responsible for the shift of the 
peak of the distribution to lower frequency. In fact we show 
below, Sec. lIVB~2l that a large population of localized states 
with n^$> n c yields a distribution which is peaked at lo w 0. 

Finally, we note that the small population of localized 
modes after the displacement also explains why the analytic 
solution Eq.d34> reproduces the exact fermionic evolution in 
Fig-El oni y qualitatively. In fact, Eq.ll34t was derived under 
the condition J2 n =n -1 Icn^ | 2 "C 1, while here n c = 12 
andE^nl^l^O.l. 



2. Localized dynamics 

In analogy to recent experiments 1 19], in this section we 
study the dipole dynamics when the depth V of the optical 
lattice is varied, while the parabolic confinement is kept con- 
stant. Then, both J and U change as a function of the lattice 
depth, as explained in Sec. I1VI The parabolic confinement 
lot = 2ir x 150Hz has been chosen to be the same as in 
Sec. llVI so that the energy spectrum exactly matches the one 
discussed there, when V is varied. 

In Fig. ^]lines and symbols correspond to the time evolu- 
tion of the interacting bosons and ideal fermions, respectively. 
In particular, the dashed-dotted, dashed and small-dotted lines 
are for bosons, while boxes, large-dots and triangles are for 
fermions with V /Er = 7, 12 and 15, respectively. For such 
lattice depths 7 = 14, 50 and 100, respectively, and the system 
is fermionized, as discussed in Sec. IIV Al As expected, the 
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FIG. 14: Center of mass position (in lattice sites) as a function of time 
calculated for different lattice depths and a fixed lot — 2tt x 150 
Hz. The dash-dotted, dashed and small-dotted lines are the exact 
solutions for interacting bosons (e in the legend) and the boxes, large- 
dots and triangles are the solutions for ideal fermions (F in the legend 
) for Vo/Er = 7, 12 and 15, respectively. T a = 2-k/lo*. 



consequence of the large population of single-particle states 
which are localized in character. For the case V a = 7Er 
(QN/J w 0.58, 7 = 13.6 , q = 134.2 and n c w 8) before 
the displacement the system is fermionized but non-localized. 
On the other hand, after the displacement about 20 percent 
of the atoms occupy localized modes of the undisplaced po- 
tential. The population of the localized modes with n ^> n c 
can be directly linked to the presence of low-frequency peaks 
(pjik ~ 0) in the distribution of frequencies, Fig.^] Because 
80 percent of the atoms occupy non-localized modes the cen- 
ter of mass position can still relax to zero as shown in Figll4l 
For the cases V = YIEr and \bEn the Fermi energy is larger 
than E ne , and even before the displacement most states are lo- 
calized. After the displacement has taken place, about 60 and 
90 percent of the atoms occupy localized modes respectively, 
and the dynamics is overdamped. This is mirrored by the ap- 
pearance of a large peak at u>ik ~ in the probability distri- 
bution, Fig^] and by the fact that the center of mass position 
does not relax to zero as shown in Figfl4l 



V. CONCLUSIONS 
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FIG. 15: Probability distribution of frequencies for V /Er = 7, 12 
and 15 (see text). The frequencies u>ik are given in units of uj* . Be- 
cause lo* oc \f~J, ui* decreases with increasing lattice depth. 



agreement between the bosonic and fermionic solutions im- 
proves for larger 7-values, and is almost perfect for 7 = 100. 

Notably, for all the 7-values, no complete oscillation 
are observed, as the amplitudes of oscillations are strongly 
damped at very early times. The inhibition of the trans- 
port properties in the experiment here envisioned is a direct 



We studied the spectrum and dipolar motion of interacting 
and non-interacting one-dimensional atomic gases trapped in 
an optical lattice plus a parabolic potential using the tight- 
binding approximation. We showed that the single-atom tight- 
binding Schrodinger equation can be exactly solved by map- 
ping it onto the recurrence relation satisfied by the Fourier 
coefficients of some periodic Mathieu functions. We used 
asymptotic expansions of such functions to fully characterize 
the eigenenergies and eigenmodes of the system. Our ana- 
lytic approach is complementary to previous numerical and 
semiclassical analysis for single-atom systems. The advan- 
tage is that we could explicitly calculate the corrections to the 
harmonic oscillator spectrum introduced by the lattice. The 
knowledge of these corrections allowed us also to provide an- 
alytic expressions for the modulations of the center of mass 
motion induced by the periodic potential when trapped ideal 
bosonic and fermionic gases are suddenly displaced from the 
trap center. 

By means of numerical diagonalizations of the Bose- 
Hubbard Hamiltonian we studied the interacting many-body 
bosonic problem. First, we characterized the changes in the 
low-energy excitation spectrum as a function of lattice depth, 
by comparing it with the ideal Bose and Fermi spectra. Then, 
we stated the necessary conditions for fermionization to occur 
and showed that it takes place for a large range of experimen- 
tally accessible parameters. We clarified the required condi- 
tions for the formation of a Mott insulator at the trap center 
and linked its appearance to the population of localized states 
at the Fermi level of the correspondent ideal fermionic system. 
We then studied the many-body dipole dynamics and showed 
that, in the parameter regime where the system is expected 
to be fermionized, the knowledge of the single-particle solu- 
tions is a powerful tool for the understanding of the strongly 
correlated dynamics. By studying the distribution of the fre- 
quencies pertaining the many-body modes excited during the 
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dipole dynamics, we explicitly showed the connection be- 
tween the population of single-particle localized states with 
the inhibition of the transport properties of the system. These 
spectral analysis allowed us also to gain some insight into the 
dynamics in the weakly-interacting regime, where an analysis 
beyond mean-field is required, and in the complex interme- 
diate regime, where no mapping to single-particle solution is 
possible. 
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Note: We note that some of the points discussed here have 
been just recently pointed out in Ref.|43] where the authors 
discuss the dipole oscillations of ID bosons in the hard-core 
limit. 
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